Fit temperature models and predict growing season temperature

Author

Max Lindmark, Jan Ohlberger, Anna Gårdmark

Published

June 30, 2023

Load libraries

Code
pkgs <- c("here","tidyverse", "tidylog", "RColorBrewer", "viridis", "sdmTMB", "sdmTMBextra", "patchwork", "RCurl", "tidylog") 

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))
#> here() starts at /Users/maxlindmark/Dropbox/Max work/R/perch-growth
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr     1.1.2     ✔ readr     2.1.4
#> ✔ forcats   1.0.0     ✔ stringr   1.5.0
#> ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
#> ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
#> ✔ purrr     1.0.1     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#> 
#> Attaching package: 'tidylog'
#> 
#> 
#> The following objects are masked from 'package:dplyr':
#> 
#>     add_count, add_tally, anti_join, count, distinct, distinct_all,
#>     distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
#>     full_join, group_by, group_by_all, group_by_at, group_by_if,
#>     inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
#>     relocate, rename, rename_all, rename_at, rename_if, rename_with,
#>     right_join, sample_frac, sample_n, select, select_all, select_at,
#>     select_if, semi_join, slice, slice_head, slice_max, slice_min,
#>     slice_sample, slice_tail, summarise, summarise_all, summarise_at,
#>     summarise_if, summarize, summarize_all, summarize_at, summarize_if,
#>     tally, top_frac, top_n, transmute, transmute_all, transmute_at,
#>     transmute_if, ungroup
#> 
#> 
#> The following objects are masked from 'package:tidyr':
#> 
#>     drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
#>     spread, uncount
#> 
#> 
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> 
#> 
#> Loading required package: viridisLite
#> 
#> 
#> Attaching package: 'sdmTMBextra'
#> 
#> 
#> The following objects are masked from 'package:sdmTMB':
#> 
#>     dharma_residuals, extract_mcmc
#> 
#> 
#> 
#> Attaching package: 'RCurl'
#> 
#> 
#> The following object is masked from 'package:tidyr':
#> 
#>     complete

# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path:
home <- here::here()

Load cache

Code
# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/01-fit-temp-models-predict_cache/html"))

Read data

Code
d <- read_csv(paste0(home, "/output/temp_data_for_fitting.csv"))
#> Rows: 98616 Columns: 13
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (4): area, source, db, id
#> dbl  (8): year, temp, yday, month, day, VkDag, stn_nr, section_nr
#> date (1): date
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

d <- d |> mutate(area = as.factor(area),
                 source_f = as.factor(source),
                 year_f = as.factor(year))
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA

# Keep track of the years for which we have cohorts for matching with cohort data
gdat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv")
#> Rows: 338460 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (6): age_ring, area, gear, ID, sex, keep
#> dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

gdat$area_cohort_age <- as.factor(paste(gdat$area, gdat$cohort, gdat$age_bc))

gdat <- gdat |>
  group_by(area_cohort_age) |> 
  filter(n() > 10)
#> group_by: one grouping variable (area_cohort_age)
#> filter (grouped): removed 5,190 rows (2%), 333,270 rows remaining

gdat <- gdat |> 
  filter(age_catch > 3) |> 
  group_by(area) |>
  summarise(min = min(cohort),
            max = max(cohort)) |> 
  arrange(min)
#> filter (grouped): removed 107,058 rows (32%), 226,212 rows remaining
#> group_by: one grouping variable (area)
#> summarise: now 12 rows and 3 columns, ungrouped

d <- left_join(d, gdat, by = "area") |>
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (     0)
#>            > matched rows     98,616
#>            >                 ========
#>            > rows total       98,616
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA

# Drop data in SI_HA and BT before onset of warming?
d2 <- d |>
  mutate(discard = "N",
         discard = ifelse(area == "BT" & year <= 1980, "Y", discard),
         discard = ifelse(area == "SI_HA" & year <= 1972, "Y", discard)) |> 
  filter(discard == "N")
#> mutate: new variable 'discard' (character) with 2 unique values and 0% NA
#> filter: removed 1,146 rows (1%), 97,470 rows remaining

Fit models

Source as independent or interactive effect

Code
m <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 5),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

Try alternative degrees of freedom not evaluated

Code
#df=20 (effectively gaussian)
m2 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 20),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# df = 9
m3 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 9),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# df=4
m4 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 4),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# Plot all residuals
mcmc_res <- residuals(m, type = "mle-mcmc",
                      mcmc_samples = sdmTMBextra::predict_mle_mcmc(m,
                                                                   mcmc_iter = 201,
                                                                   mcmc_warmup = 200))

mcmc_res20 <- residuals(m2, type = "mle-mcmc",
                        mcmc_samples = sdmTMBextra::predict_mle_mcmc(m2,
                                                                     mcmc_iter = 201,
                                                                     mcmc_warmup = 200))

mcmc_res9 <- residuals(m3, type = "mle-mcmc",
                       mcmc_samples = sdmTMBextra::predict_mle_mcmc(m3,
                                                                    mcmc_iter = 201,
                                                                    mcmc_warmup = 200))

mcmc_res4 <- residuals(m4, type = "mle-mcmc",
                       mcmc_samples = sdmTMBextra::predict_mle_mcmc(m4,
                                                                    mcmc_iter = 201,
                                                                    mcmc_warmup = 200))

dres <- d |> mutate("df=5" = mcmc_res,
                    "df=20" = mcmc_res20,
                    "df=9" = mcmc_res9,
                    "df=4" = mcmc_res4) |> 
  select(`df=5`, `df=20`, `df=9`, `df=4`) |> 
  pivot_longer(everything())

ggplot(dres, aes(sample = value)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~factor(name, levels = c("df=20", "df=9", "df=5", "df=4"))) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") + 
  theme(aspect.ratio = 1)

Check fit

Code
sanity(m)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

mcmc_res <- residuals(m, type = "mle-mcmc",
                      mcmc_samples = sdmTMBextra::predict_mle_mcmc(m,
                                                                   mcmc_iter = 201,
                                                                   mcmc_warmup = 200))
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.021513 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 215.13 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 259.439 seconds (Warm-up)
#> Chain 1:                0.256778 seconds (Sampling)
#> Chain 1:                259.696 seconds (Total)
#> Chain 1:

ggplot(d, aes(sample = mcmc_res)) +
  stat_qq() +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") + 
  theme(aspect.ratio = 1)

Code

ggsave(paste0(home, "/figures/supp/full_model/qq_temp.pdf"), width = 11, height = 11, units = "cm")

summary(m)
#> Model fit by ML ['sdmTMB']
#> Formula: temp ~ area * year_f + source_f + s(yday, by = area, bs = "cc")
#> Mesh: NULL (isotropic covariance)
#> Data: d
#> Family: student(link = 'identity')
#>  
#>                      coef.est coef.se
#> (Intercept)              8.03    0.50
#> areaBT                  -0.38    0.77
#> areaFB                  -0.04    0.78
#> areaFM                   0.05    0.82
#> areaHO                  -0.36    0.85
#> areaJM                   0.58    0.80
#> areaMU                   1.32    0.76
#> areaRA                  -1.05    0.84
#> areaSI_EK                2.03    0.78
#> areaSI_HA                1.94    0.77
#> areaTH                   2.09    0.71
#> areaVN                   1.44    0.71
#> year_f1941              -0.31    0.73
#> year_f1942              -0.30    0.71
#> year_f1943               1.36    0.72
#> year_f1944               1.36    0.72
#> year_f1945               1.69    0.71
#> year_f1946               1.93    0.74
#> year_f1947               0.87    0.73
#> year_f1948               1.41    0.71
#> year_f1949               1.38    0.71
#> year_f1950               1.16    0.70
#> year_f1951               1.05    0.72
#> year_f1952               0.22    0.70
#> year_f1953               1.63    0.71
#> year_f1954               1.61    0.71
#> year_f1955               1.57    0.71
#> year_f1956               1.05    0.71
#> year_f1957               1.25    0.71
#> year_f1958               1.12    0.72
#> year_f1959               2.26    0.71
#> year_f1960               1.28    0.69
#> year_f1961               2.32    0.72
#> year_f1962               1.38    0.73
#> year_f1963               1.59    0.71
#> year_f1964               1.30    0.72
#> year_f1965               1.23    0.71
#> year_f1966               1.08    0.70
#> year_f1967               1.92    0.70
#> year_f1968               1.62    0.70
#> year_f1969               1.27    0.69
#> year_f1970               1.25    0.69
#> year_f1971               1.19    0.69
#> year_f1972               1.42    0.71
#> year_f1973               2.01    0.74
#> year_f1974               1.86    0.69
#> year_f1975               2.22    0.70
#> year_f1976               0.61    0.70
#> year_f1977               0.61    0.70
#> year_f1978               0.47    0.69
#> year_f1979               0.10    0.71
#> year_f1980               0.57    0.72
#> year_f1981               0.54    0.70
#> year_f1982               0.74    0.71
#> year_f1983               1.30    0.70
#> year_f1984               1.29    0.70
#> year_f1985              -0.11    0.70
#> year_f1986               0.38    0.71
#> year_f1987              -0.82    0.71
#> year_f1988               1.33    0.71
#> year_f1989               2.11    0.70
#> year_f1990               2.33    0.72
#> year_f1991               0.84    0.58
#> year_f1992               0.25    0.55
#> year_f1993               0.03    0.54
#> year_f1994               0.44    0.53
#> year_f1995               1.48    0.55
#> year_f1996              -0.11    0.54
#> year_f1997               3.72    0.55
#> year_f1998              -0.23    0.51
#> year_f1999               1.40    0.52
#> year_f2000               1.12    0.51
#> year_f2001               1.87    0.51
#> year_f2002               3.10    0.52
#> year_f2003               2.44    0.52
#> year_f2004               1.25    0.51
#> year_f2005               1.93    0.51
#> year_f2006               0.11    0.57
#> year_f2007               0.58    0.56
#> year_f2008               1.45    0.56
#> year_f2009               1.77    0.55
#> year_f2010               1.33    0.73
#> year_f2011               1.75    0.73
#> year_f2012               1.74    0.70
#> year_f2013               2.14    0.72
#> year_f2014               2.50    0.71
#> year_f2015               2.55    0.70
#> year_f2016               2.65    0.70
#> year_f2017               2.20    0.69
#> year_f2018               2.80    0.71
#> year_f2019               2.14    0.70
#> year_f2020               3.22    0.69
#> year_f2021               2.33    0.71
#> year_f2022               2.11    1.29
#> source_ffishing          1.55    0.05
#> source_flogger           1.26    0.04
#> areaBT:year_f1941       -0.09    1.15
#> areaFB:year_f1941       -0.01    1.18
#> areaFM:year_f1941       -0.01    1.25
#> areaHO:year_f1941        0.41    1.25
#> areaJM:year_f1941       -0.28    1.21
#> areaMU:year_f1941       -0.15    1.13
#> areaRA:year_f1941        0.11    1.22
#> areaSI_EK:year_f1941    -0.43    1.16
#> areaSI_HA:year_f1941    -0.45    1.16
#> areaTH:year_f1941       -0.56    1.05
#> areaVN:year_f1941       -0.46    1.03
#> areaBT:year_f1942       -0.02    1.10
#> areaFB:year_f1942        0.14    1.11
#> areaFM:year_f1942        0.15    1.17
#> areaHO:year_f1942        0.38    1.19
#> areaJM:year_f1942       -0.03    1.17
#> areaMU:year_f1942       -0.17    1.07
#> areaRA:year_f1942        0.08    1.17
#> areaSI_EK:year_f1942    -0.18    1.13
#> areaSI_HA:year_f1942    -0.17    1.16
#> areaTH:year_f1942        0.08    1.05
#> areaVN:year_f1942       -0.25    1.04
#> areaBT:year_f1943       -0.09    1.10
#> areaFB:year_f1943       -0.53    1.12
#> areaFM:year_f1943       -0.72    1.17
#> areaHO:year_f1943       -0.81    1.20
#> areaJM:year_f1943       -0.59    1.15
#> areaMU:year_f1943       -0.40    1.07
#> areaRA:year_f1943       -0.72    1.21
#> areaSI_EK:year_f1943    -0.31    1.09
#> areaSI_HA:year_f1943    -0.34    1.10
#> areaTH:year_f1943        0.03    1.02
#> areaVN:year_f1943       -0.03    1.02
#> areaBT:year_f1944        0.17    1.14
#> areaFB:year_f1944       -0.31    1.15
#> areaFM:year_f1944       -0.32    1.23
#> areaHO:year_f1944       -0.40    1.24
#> areaJM:year_f1944       -0.43    1.18
#> areaMU:year_f1944       -0.25    1.08
#> areaRA:year_f1944       -0.67    1.22
#> areaSI_EK:year_f1944    -0.17    1.11
#> areaSI_HA:year_f1944    -0.22    1.12
#> areaTH:year_f1944       -0.20    1.02
#> areaVN:year_f1944       -0.20    1.02
#> areaBT:year_f1945        0.18    1.07
#> areaFB:year_f1945       -0.17    1.09
#> areaFM:year_f1945       -0.35    1.15
#> areaHO:year_f1945       -0.68    1.18
#> areaJM:year_f1945       -0.14    1.13
#> areaMU:year_f1945       -0.16    1.05
#> areaRA:year_f1945       -1.07    1.19
#> areaSI_EK:year_f1945     0.07    1.08
#> areaSI_HA:year_f1945     0.03    1.08
#> areaTH:year_f1945        0.27    1.00
#> areaVN:year_f1945        0.24    1.00
#> areaBT:year_f1946        0.32    1.12
#> areaFB:year_f1946        0.20    1.11
#> areaFM:year_f1946       -0.01    1.17
#> areaHO:year_f1946       -0.55    1.20
#> areaJM:year_f1946        0.24    1.15
#> areaMU:year_f1946       -0.08    1.06
#> areaRA:year_f1946       -1.61    1.19
#> areaSI_EK:year_f1946    -0.01    1.10
#> areaSI_HA:year_f1946    -0.12    1.11
#> areaTH:year_f1946       -0.07    1.03
#> areaVN:year_f1946        0.16    1.04
#> areaBT:year_f1947        0.30    1.11
#> areaFB:year_f1947       -0.03    1.11
#> areaFM:year_f1947       -0.16    1.17
#> areaHO:year_f1947       -0.95    1.19
#> areaJM:year_f1947        0.11    1.16
#> areaMU:year_f1947       -0.01    1.07
#> areaRA:year_f1947       -1.79    1.22
#> areaSI_EK:year_f1947     0.20    1.13
#> areaSI_HA:year_f1947     0.22    1.15
#> areaTH:year_f1947        0.21    1.06
#> areaVN:year_f1947        0.36    1.05
#> areaBT:year_f1948        0.03    1.09
#> areaFB:year_f1948       -0.04    1.10
#> areaFM:year_f1948       -0.09    1.16
#> areaHO:year_f1948       -0.19    1.19
#> areaJM:year_f1948       -0.08    1.13
#> areaMU:year_f1948       -0.26    1.06
#> areaRA:year_f1948       -0.76    1.17
#> areaSI_EK:year_f1948    -0.24    1.10
#> areaSI_HA:year_f1948    -0.22    1.10
#> areaTH:year_f1948       -0.27    1.01
#> areaVN:year_f1948       -0.01    1.00
#> areaBT:year_f1949        0.10    1.10
#> areaFB:year_f1949       -0.18    1.14
#> areaFM:year_f1949       -0.34    1.20
#> areaHO:year_f1949       -0.55    1.25
#> areaJM:year_f1949       -0.18    1.17
#> areaMU:year_f1949        0.06    1.09
#> areaRA:year_f1949       -1.36    1.23
#> areaSI_EK:year_f1949     0.31    1.14
#> areaSI_HA:year_f1949     0.23    1.14
#> areaTH:year_f1949        0.44    1.04
#> areaVN:year_f1949        0.43    1.01
#> areaBT:year_f1950        0.10    1.08
#> areaFB:year_f1950        0.03    1.11
#> areaFM:year_f1950       -0.11    1.18
#> areaHO:year_f1950       -0.22    1.20
#> areaJM:year_f1950        0.23    1.14
#> areaMU:year_f1950        0.02    1.05
#> areaRA:year_f1950       -0.58    1.22
#> areaSI_EK:year_f1950     0.24    1.08
#> areaSI_HA:year_f1950     0.16    1.08
#> areaTH:year_f1950        0.28    1.00
#> areaVN:year_f1950        0.39    1.00
#> areaBT:year_f1951        0.10    1.14
#> areaFB:year_f1951       -0.29    1.17
#> areaFM:year_f1951       -0.42    1.25
#> areaHO:year_f1951       -0.45    1.27
#> areaJM:year_f1951       -0.34    1.21
#> areaMU:year_f1951       -0.03    1.11
#> areaRA:year_f1951       -1.15    1.24
#> areaSI_EK:year_f1951     0.11    1.16
#> areaSI_HA:year_f1951     0.06    1.17
#> areaTH:year_f1951        0.00    1.05
#> areaVN:year_f1951        0.14    1.02
#> areaBT:year_f1952       -0.05    1.06
#> areaFB:year_f1952       -0.07    1.08
#> areaFM:year_f1952       -0.18    1.13
#> areaHO:year_f1952       -0.39    1.19
#> areaJM:year_f1952        0.18    1.11
#> areaMU:year_f1952       -0.08    1.06
#> areaRA:year_f1952       -0.91    1.17
#> areaSI_EK:year_f1952     0.10    1.09
#> areaSI_HA:year_f1952     0.08    1.08
#> areaTH:year_f1952        0.33    1.01
#> areaVN:year_f1952        0.40    1.01
#> areaBT:year_f1953        0.23    1.06
#> areaFB:year_f1953        0.17    1.10
#> areaFM:year_f1953        0.09    1.14
#> areaHO:year_f1953       -0.30    1.20
#> areaJM:year_f1953        0.01    1.14
#> areaMU:year_f1953       -0.12    1.06
#> areaRA:year_f1953       -1.12    1.18
#> areaSI_EK:year_f1953    -0.30    1.09
#> areaSI_HA:year_f1953    -0.33    1.09
#> areaTH:year_f1953       -0.17    1.01
#> areaVN:year_f1953       -0.05    1.01
#> areaBT:year_f1954        0.17    1.05
#> areaFB:year_f1954        0.20    1.10
#> areaFM:year_f1954        0.13    1.15
#> areaHO:year_f1954       -0.15    1.20
#> areaJM:year_f1954       -0.09    1.13
#> areaMU:year_f1954       -0.37    1.09
#> areaRA:year_f1954       -1.28    1.16
#> areaSI_EK:year_f1954    -0.64    1.11
#> areaSI_HA:year_f1954    -0.60    1.10
#> areaTH:year_f1954       -0.59    1.02
#> areaVN:year_f1954       -0.38    1.00
#> areaBT:year_f1955        0.31    1.09
#> areaFB:year_f1955        0.63    1.13
#> areaFM:year_f1955        0.74    1.17
#> areaHO:year_f1955        0.23    1.22
#> areaJM:year_f1955        0.50    1.18
#> areaMU:year_f1955        0.08    1.10
#> areaRA:year_f1955       -1.47    1.17
#> areaSI_EK:year_f1955     0.03    1.14
#> areaSI_HA:year_f1955    -0.06    1.15
#> areaTH:year_f1955       -0.28    1.06
#> areaVN:year_f1955       -0.34    1.04
#> areaBT:year_f1956        0.05    1.06
#> areaFB:year_f1956        0.08    1.11
#> areaFM:year_f1956       -0.03    1.15
#> areaHO:year_f1956       -0.04    1.23
#> areaJM:year_f1956       -0.28    1.13
#> areaMU:year_f1956       -0.60    1.09
#> areaRA:year_f1956       -0.99    1.17
#> areaSI_EK:year_f1956    -0.90    1.11
#> areaSI_HA:year_f1956    -0.86    1.10
#> areaTH:year_f1956       -0.76    1.02
#> areaVN:year_f1956       -0.65    1.01
#> areaBT:year_f1957        0.18    1.06
#> areaFB:year_f1957        0.40    1.12
#> areaFM:year_f1957        0.40    1.16
#> areaHO:year_f1957        0.30    1.20
#> areaJM:year_f1957        0.26    1.14
#> areaMU:year_f1957       -0.25    1.11
#> areaRA:year_f1957       -0.86    1.18
#> areaSI_EK:year_f1957    -0.33    1.11
#> areaSI_HA:year_f1957    -0.30    1.10
#> areaTH:year_f1957       -0.11    1.02
#> areaVN:year_f1957       -0.08    1.01
#> areaBT:year_f1958        0.16    1.08
#> areaFB:year_f1958        0.36    1.18
#> areaFM:year_f1958        0.33    1.23
#> areaHO:year_f1958        0.44    1.25
#> areaJM:year_f1958        0.03    1.22
#> areaMU:year_f1958       -0.15    1.16
#> areaRA:year_f1958       -0.79    1.21
#> areaSI_EK:year_f1958    -0.26    1.20
#> areaSI_HA:year_f1958    -0.23    1.18
#> areaTH:year_f1958       -0.17    1.07
#> areaVN:year_f1958       -0.20    1.02
#> areaBT:year_f1959        0.17    1.05
#> areaFB:year_f1959        0.32    1.08
#> areaFM:year_f1959        0.24    1.13
#> areaHO:year_f1959       -0.14    1.17
#> areaJM:year_f1959        0.25    1.12
#> areaMU:year_f1959       -0.17    1.07
#> areaRA:year_f1959       -1.56    1.16
#> areaSI_EK:year_f1959    -0.16    1.09
#> areaSI_HA:year_f1959    -0.23    1.09
#> areaTH:year_f1959       -0.06    1.01
#> areaVN:year_f1959       -0.04    1.01
#> areaBT:year_f1960        0.24    1.05
#> areaFB:year_f1960        0.26    1.08
#> areaFM:year_f1960        0.23    1.13
#> areaHO:year_f1960       -0.11    1.17
#> areaJM:year_f1960        0.11    1.12
#> areaMU:year_f1960       -0.18    1.06
#> areaRA:year_f1960       -0.96    1.15
#> areaSI_EK:year_f1960    -0.27    1.09
#> areaSI_HA:year_f1960    -0.26    1.09
#> areaTH:year_f1960       -0.15    1.00
#> areaVN:year_f1960       -0.08    0.99
#> areaBT:year_f1961        0.10    1.06
#> areaFB:year_f1961        0.04    1.14
#> areaFM:year_f1961       -0.13    1.18
#> areaHO:year_f1961       -0.57    1.22
#> areaJM:year_f1961       -0.15    1.17
#> areaMU:year_f1961       -0.20    1.12
#> areaRA:year_f1961       -1.59    1.20
#> areaSI_EK:year_f1961    -0.46    1.13
#> areaSI_HA:year_f1961    -0.49    1.11
#> areaTH:year_f1961       -0.40    1.03
#> areaVN:year_f1961       -0.22    1.01
#> areaBT:year_f1962        0.02    1.07
#> areaFB:year_f1962        0.17    1.17
#> areaFM:year_f1962       -0.02    1.22
#> areaHO:year_f1962       -0.08    1.26
#> areaJM:year_f1962       -0.18    1.20
#> areaMU:year_f1962       -0.41    1.17
#> areaRA:year_f1962       -1.14    1.21
#> areaSI_EK:year_f1962    -0.88    1.19
#> areaSI_HA:year_f1962    -0.81    1.17
#> areaTH:year_f1962       -0.80    1.07
#> areaVN:year_f1962       -0.51    1.03
#> areaBT:year_f1963        0.22    1.08
#> areaFB:year_f1963        0.13    1.10
#> areaFM:year_f1963        0.00    1.16
#> areaHO:year_f1963       -0.23    1.20
#> areaJM:year_f1963       -0.22    1.15
#> areaMU:year_f1963       -0.30    1.06
#> areaRA:year_f1963       -1.02    1.19
#> areaSI_EK:year_f1963    -0.65    1.10
#> areaSI_HA:year_f1963    -0.68    1.12
#> areaTH:year_f1963       -0.89    1.03
#> areaVN:year_f1963       -0.54    1.01
#> areaBT:year_f1964        0.04    1.06
#> areaFB:year_f1964        0.07    1.13
#> areaFM:year_f1964       -0.13    1.17
#> areaHO:year_f1964       -0.05    1.22
#> areaJM:year_f1964       -0.27    1.15
#> areaMU:year_f1964       -0.37    1.11
#> areaRA:year_f1964       -0.89    1.21
#> areaSI_EK:year_f1964    -0.67    1.13
#> areaSI_HA:year_f1964    -0.64    1.12
#> areaTH:year_f1964       -0.31    1.03
#> areaVN:year_f1964       -0.33    1.01
#> areaBT:year_f1965       -0.01    1.06
#> areaFB:year_f1965        0.19    1.14
#> areaFM:year_f1965        0.03    1.19
#> areaHO:year_f1965        0.19    1.23
#> areaJM:year_f1965       -0.05    1.17
#> areaMU:year_f1965       -0.46    1.12
#> areaRA:year_f1965       -0.84    1.18
#> areaSI_EK:year_f1965    -0.82    1.14
#> areaSI_HA:year_f1965    -0.84    1.13
#> areaTH:year_f1965       -0.67    1.04
#> areaVN:year_f1965       -0.58    1.01
#> areaBT:year_f1966        0.20    1.07
#> areaFB:year_f1966        0.12    1.09
#> areaFM:year_f1966        0.08    1.13
#> areaHO:year_f1966       -0.33    1.18
#> areaJM:year_f1966       -0.02    1.12
#> areaMU:year_f1966       -0.10    1.06
#> areaRA:year_f1966       -1.18    1.15
#> areaSI_EK:year_f1966    -0.18    1.09
#> areaSI_HA:year_f1966    -0.15    1.09
#> areaTH:year_f1966       -0.09    1.02
#> areaVN:year_f1966       -0.06    1.01
#> areaBT:year_f1967        0.15    1.06
#> areaFB:year_f1967        0.22    1.12
#> areaFM:year_f1967        0.12    1.17
#> areaHO:year_f1967        0.01    1.21
#> areaJM:year_f1967        0.02    1.17
#> areaMU:year_f1967       -0.11    1.08
#> areaRA:year_f1967       -1.02    1.20
#> areaSI_EK:year_f1967    -0.23    1.11
#> areaSI_HA:year_f1967    -0.29    1.11
#> areaTH:year_f1967       -0.15    1.01
#> areaVN:year_f1967       -0.10    1.00
#> areaBT:year_f1968        0.11    1.04
#> areaFB:year_f1968        0.25    1.10
#> areaFM:year_f1968        0.17    1.14
#> areaHO:year_f1968        0.11    1.20
#> areaJM:year_f1968        0.10    1.14
#> areaMU:year_f1968       -0.20    1.06
#> areaRA:year_f1968       -1.16    1.19
#> areaSI_EK:year_f1968    -0.32    1.09
#> areaSI_HA:year_f1968    -0.39    1.09
#> areaTH:year_f1968       -0.18    1.01
#> areaVN:year_f1968       -0.17    0.99
#> areaBT:year_f1969        0.30    1.06
#> areaFB:year_f1969        0.39    1.08
#> areaFM:year_f1969        0.37    1.14
#> areaHO:year_f1969        0.03    1.18
#> areaJM:year_f1969        0.39    1.14
#> areaMU:year_f1969        0.06    1.05
#> areaRA:year_f1969       -0.96    1.16
#> areaSI_EK:year_f1969     0.23    1.09
#> areaSI_HA:year_f1969     0.19    1.10
#> areaTH:year_f1969        0.23    1.01
#> areaVN:year_f1969        0.16    1.00
#> areaBT:year_f1970        0.19    1.05
#> areaFB:year_f1970        0.13    1.08
#> areaFM:year_f1970        0.05    1.13
#> areaHO:year_f1970       -0.09    1.18
#> areaJM:year_f1970       -0.02    1.12
#> areaMU:year_f1970       -0.21    1.05
#> areaRA:year_f1970       -0.81    1.16
#> areaSI_EK:year_f1970    -0.28    1.08
#> areaSI_HA:year_f1970    -0.28    1.09
#> areaTH:year_f1970       -0.27    1.01
#> areaVN:year_f1970       -0.16    0.99
#> areaBT:year_f1971        0.24    1.05
#> areaFB:year_f1971        0.23    1.10
#> areaFM:year_f1971        0.20    1.15
#> areaHO:year_f1971       -0.25    1.19
#> areaJM:year_f1971        0.27    1.13
#> areaMU:year_f1971       -0.02    1.08
#> areaRA:year_f1971       -1.19    1.14
#> areaSI_EK:year_f1971     0.07    1.10
#> areaSI_HA:year_f1971     0.08    1.10
#> areaTH:year_f1971        0.19    1.01
#> areaVN:year_f1971        0.22    0.99
#> areaBT:year_f1972        0.31    1.11
#> areaFB:year_f1972        0.52    1.11
#> areaFM:year_f1972        0.53    1.18
#> areaHO:year_f1972       -0.28    1.19
#> areaJM:year_f1972       -0.40    1.01
#> areaMU:year_f1972        0.40    1.07
#> areaRA:year_f1972       -0.91    1.17
#> areaSI_EK:year_f1972     0.42    1.10
#> areaSI_HA:year_f1972     0.99    0.95
#> areaTH:year_f1972        0.26    1.01
#> areaVN:year_f1972        0.39    1.02
#> areaBT:year_f1973        0.07    1.08
#> areaFB:year_f1973        0.29    1.09
#> areaFM:year_f1973        0.25    1.13
#> areaHO:year_f1973        0.07    1.18
#> areaJM:year_f1973       -3.61    1.03
#> areaMU:year_f1973       -0.19    1.10
#> areaRA:year_f1973       -0.81    1.15
#> areaSI_EK:year_f1973    -0.11    1.10
#> areaSI_HA:year_f1973     2.15    0.97
#> areaTH:year_f1973       -0.01    1.04
#> areaVN:year_f1973        0.14    1.05
#> areaBT:year_f1974        0.30    1.05
#> areaFB:year_f1974        0.32    1.10
#> areaFM:year_f1974        0.21    1.14
#> areaHO:year_f1974        0.12    1.20
#> areaJM:year_f1974       -4.66    0.98
#> areaMU:year_f1974       -0.06    1.07
#> areaRA:year_f1974       -1.19    1.16
#> areaSI_EK:year_f1974    -0.29    1.09
#> areaSI_HA:year_f1974    -0.58    0.93
#> areaTH:year_f1974       -0.39    1.00
#> areaVN:year_f1974       -0.11    0.99
#> areaBT:year_f1975        0.20    1.05
#> areaFB:year_f1975        0.48    1.10
#> areaFM:year_f1975       -1.39    1.04
#> areaHO:year_f1975        0.49    1.19
#> areaJM:year_f1975       -2.72    0.99
#> areaMU:year_f1975        0.02    1.07
#> areaRA:year_f1975       -0.60    1.16
#> areaSI_EK:year_f1975     0.10    1.10
#> areaSI_HA:year_f1975     5.64    0.93
#> areaTH:year_f1975        0.02    1.02
#> areaVN:year_f1975        0.09    1.00
#> areaBT:year_f1976        0.24    1.05
#> areaFB:year_f1976       -1.06    1.01
#> areaFM:year_f1976       -0.01    1.03
#> areaHO:year_f1976        0.51    1.18
#> areaJM:year_f1976       -0.77    0.99
#> areaMU:year_f1976        0.05    1.07
#> areaRA:year_f1976       -0.43    1.17
#> areaSI_EK:year_f1976     0.36    1.09
#> areaSI_HA:year_f1976     6.18    0.93
#> areaTH:year_f1976        0.63    1.00
#> areaVN:year_f1976        0.39    1.00
#> areaBT:year_f1977        0.21    1.05
#> areaFB:year_f1977       -2.77    0.99
#> areaFM:year_f1977       -1.61    1.01
#> areaHO:year_f1977       -0.09    1.21
#> areaJM:year_f1977       -2.36    0.98
#> areaMU:year_f1977        0.08    1.08
#> areaRA:year_f1977       -0.72    1.20
#> areaSI_EK:year_f1977     0.10    1.09
#> areaSI_HA:year_f1977     5.55    0.93
#> areaTH:year_f1977        0.32    1.01
#> areaVN:year_f1977        0.27    0.99
#> areaBT:year_f1978        0.28    1.05
#> areaFB:year_f1978       -2.08    0.97
#> areaFM:year_f1978       -1.50    1.05
#> areaHO:year_f1978        0.00    1.19
#> areaJM:year_f1978       -1.90    0.99
#> areaMU:year_f1978        0.01    1.08
#> areaRA:year_f1978       -0.91    1.14
#> areaSI_EK:year_f1978     0.01    1.10
#> areaSI_HA:year_f1978     5.74    0.93
#> areaTH:year_f1978        0.27    1.01
#> areaVN:year_f1978        0.22    0.99
#> areaBT:year_f1979        0.18    1.10
#> areaFB:year_f1979       -0.91    0.99
#> areaFM:year_f1979       -0.88    1.05
#> areaHO:year_f1979       -0.14    1.19
#> areaJM:year_f1979       -3.64    1.00
#> areaMU:year_f1979        0.13    1.05
#> areaRA:year_f1979       -0.54    1.16
#> areaSI_EK:year_f1979     0.25    1.09
#> areaSI_HA:year_f1979     4.12    0.94
#> areaTH:year_f1979        0.24    1.02
#> areaVN:year_f1979        0.34    1.01
#> areaBT:year_f1980        0.26    1.11
#> areaFB:year_f1980       -1.54    1.04
#> areaFM:year_f1980       -0.18    1.07
#> areaHO:year_f1980        0.12    1.22
#> areaJM:year_f1980       -1.07    0.99
#> areaMU:year_f1980        0.03    1.06
#> areaRA:year_f1980       -0.58    1.18
#> areaSI_EK:year_f1980     0.12    1.09
#> areaSI_HA:year_f1980     5.85    0.95
#> areaTH:year_f1980        0.10    1.01
#> areaVN:year_f1980        0.18    1.01
#> areaBT:year_f1981        0.23    1.06
#> areaFB:year_f1981       -1.60    1.00
#> areaFM:year_f1981       -0.36    1.05
#> areaHO:year_f1981       -0.35    1.20
#> areaJM:year_f1981       -1.48    0.98
#> areaMU:year_f1981        0.14    1.06
#> areaRA:year_f1981       -1.26    1.17
#> areaSI_EK:year_f1981     0.16    1.09
#> areaSI_HA:year_f1981     4.75    0.93
#> areaTH:year_f1981        0.23    1.00
#> areaVN:year_f1981        0.28    0.99
#> areaBT:year_f1982        0.12    1.11
#> areaFB:year_f1982       -2.22    1.00
#> areaFM:year_f1982       -1.08    1.04
#> areaHO:year_f1982       -0.51    1.23
#> areaJM:year_f1982       -3.30    0.99
#> areaMU:year_f1982        0.11    1.07
#> areaRA:year_f1982       -0.63    1.20
#> areaSI_EK:year_f1982     0.25    1.11
#> areaSI_HA:year_f1982     5.25    0.94
#> areaTH:year_f1982        0.49    1.03
#> areaVN:year_f1982        0.45    1.01
#> areaBT:year_f1983        0.31    1.04
#> areaFB:year_f1983       -2.15    0.97
#> areaFM:year_f1983       -0.80    0.98
#> areaHO:year_f1983       -0.13    1.17
#> areaJM:year_f1983       -3.19    0.98
#> areaMU:year_f1983        0.02    1.06
#> areaRA:year_f1983       -0.89    1.14
#> areaSI_EK:year_f1983     0.10    1.08
#> areaSI_HA:year_f1983     5.26    0.93
#> areaTH:year_f1983        0.17    1.00
#> areaVN:year_f1983        0.21    1.00
#> areaBT:year_f1984        6.63    0.93
#> areaFB:year_f1984       -0.52    1.01
#> areaFM:year_f1984        0.40    0.98
#> areaHO:year_f1984       -0.34    1.18
#> areaJM:year_f1984       -1.47    0.98
#> areaMU:year_f1984        0.02    1.06
#> areaRA:year_f1984       -1.18    1.16
#> areaSI_EK:year_f1984     0.02    1.10
#> areaSI_HA:year_f1984     6.28    0.94
#> areaTH:year_f1984       -0.07    1.02
#> areaVN:year_f1984        0.10    0.99
#> areaBT:year_f1985        0.22    1.10
#> areaFB:year_f1985       -0.21    0.96
#> areaFM:year_f1985        2.95    0.98
#> areaHO:year_f1985        0.05    1.19
#> areaJM:year_f1985       -5.33    0.98
#> areaMU:year_f1985        0.10    1.06
#> areaRA:year_f1985       -0.42    1.19
#> areaSI_EK:year_f1985     0.29    1.09
#> areaSI_HA:year_f1985     2.87    0.95
#> areaTH:year_f1985        0.21    1.01
#> areaVN:year_f1985        0.32    1.00
#> areaBT:year_f1986        8.00    0.93
#> areaFB:year_f1986       -2.10    0.98
#> areaFM:year_f1986        8.44    0.98
#> areaHO:year_f1986       -0.16    1.16
#> areaJM:year_f1986       -4.08    1.04
#> areaMU:year_f1986       -0.06    1.06
#> areaRA:year_f1986       -0.14    1.14
#> areaSI_EK:year_f1986    -0.15    1.09
#> areaSI_HA:year_f1986     4.12    0.95
#> areaTH:year_f1986       -0.15    1.02
#> areaVN:year_f1986        0.04    1.01
#> areaBT:year_f1987        8.33    0.93
#> areaFB:year_f1987       -1.80    0.98
#> areaFM:year_f1987        4.56    1.00
#> areaHO:year_f1987       -0.21    1.25
#> areaJM:year_f1987       -1.27    0.97
#> areaMU:year_f1987        0.07    1.09
#> areaRA:year_f1987       -0.12    1.23
#> areaSI_EK:year_f1987    -2.46    0.98
#> areaSI_HA:year_f1987     5.89    0.95
#> areaTH:year_f1987        0.54    1.04
#> areaVN:year_f1987        0.55    1.01
#> areaBT:year_f1988        6.94    0.93
#> areaFB:year_f1988       -3.00    0.97
#> areaFM:year_f1988        5.64    0.99
#> areaHO:year_f1988        0.18    1.17
#> areaJM:year_f1988       -3.34    0.97
#> areaMU:year_f1988       -0.06    1.07
#> areaRA:year_f1988       -1.13    1.16
#> areaSI_EK:year_f1988    -4.88    1.00
#> areaSI_HA:year_f1988     3.81    0.94
#> areaTH:year_f1988       -0.08    1.01
#> areaVN:year_f1988       -0.07    1.02
#> areaBT:year_f1989        6.64    0.92
#> areaFB:year_f1989       -4.02    0.96
#> areaFM:year_f1989       -0.79    0.97
#> areaHO:year_f1989       -2.60    1.00
#> areaJM:year_f1989       -2.40    0.95
#> areaMU:year_f1989       -0.18    1.05
#> areaRA:year_f1989       -1.08    1.15
#> areaSI_EK:year_f1989    -2.39    0.94
#> areaSI_HA:year_f1989     4.63    0.93
#> areaTH:year_f1989       -0.01    1.00
#> areaVN:year_f1989        0.03    1.00
#> areaBT:year_f1990        6.02    0.94
#> areaFB:year_f1990       -3.09    0.96
#> areaFM:year_f1990        1.86    0.98
#> areaHO:year_f1990       -2.14    1.00
#> areaJM:year_f1990        0.35    0.97
#> areaMU:year_f1990       -0.24    1.07
#> areaRA:year_f1990       -0.96    1.15
#> areaSI_EK:year_f1990    -1.97    0.94
#> areaSI_HA:year_f1990     5.25    0.94
#> areaTH:year_f1990       -0.07    1.02
#> areaVN:year_f1990        0.04    1.02
#> areaBT:year_f1991        7.53    0.83
#> areaFB:year_f1991        0.41    0.84
#> areaFM:year_f1991        3.45    0.89
#> areaHO:year_f1991       -0.57    0.90
#> areaJM:year_f1991       -0.43    0.87
#> areaMU:year_f1991       -0.28    0.86
#> areaRA:year_f1991        0.50    0.91
#> areaSI_EK:year_f1991    -0.66    0.91
#> areaSI_HA:year_f1991     5.71    0.84
#> areaTH:year_f1991        0.57    0.93
#> areaVN:year_f1991        0.74    0.92
#> areaBT:year_f1992        8.72    0.81
#> areaFB:year_f1992        0.89    0.82
#> areaFM:year_f1992        1.27    0.86
#> areaHO:year_f1992       -0.03    0.89
#> areaJM:year_f1992        1.56    0.89
#> areaMU:year_f1992        1.16    0.87
#> areaRA:year_f1992        1.92    0.92
#> areaSI_EK:year_f1992    -1.24    0.82
#> areaSI_HA:year_f1992     5.05    0.82
#> areaTH:year_f1992        1.48    0.91
#> areaVN:year_f1992        1.64    0.92
#> areaBT:year_f1993        8.84    0.80
#> areaFB:year_f1993        0.58    0.81
#> areaFM:year_f1993        0.68    0.85
#> areaHO:year_f1993       -0.88    0.88
#> areaJM:year_f1993        0.82    0.83
#> areaMU:year_f1993        0.06    0.87
#> areaRA:year_f1993        0.30    1.06
#> areaSI_EK:year_f1993    -0.60    0.82
#> areaSI_HA:year_f1993     5.64    0.82
#> areaTH:year_f1993        0.85    0.90
#> areaVN:year_f1993        0.98    0.90
#> areaBT:year_f1994        8.98    0.80
#> areaFB:year_f1994        0.93    0.81
#> areaFM:year_f1994        2.53    0.85
#> areaHO:year_f1994        1.05    0.88
#> areaJM:year_f1994        1.44    0.83
#> areaMU:year_f1994       -1.18    0.90
#> areaRA:year_f1994        0.28    0.90
#> areaSI_EK:year_f1994     0.51    0.82
#> areaSI_HA:year_f1994     6.64    0.81
#> areaTH:year_f1994        0.97    0.90
#> areaVN:year_f1994        0.93    0.90
#> areaBT:year_f1995        6.66    0.81
#> areaFB:year_f1995       -0.24    0.82
#> areaFM:year_f1995        0.75    0.86
#> areaHO:year_f1995       -1.70    0.89
#> areaJM:year_f1995        0.87    0.84
#> areaMU:year_f1995       -0.23    0.87
#> areaRA:year_f1995       -4.59    0.92
#> areaSI_EK:year_f1995    -3.68    0.83
#> areaSI_HA:year_f1995     3.96    0.82
#> areaTH:year_f1995        0.28    0.91
#> areaVN:year_f1995        0.27    0.80
#> areaBT:year_f1996        7.47    0.81
#> areaFB:year_f1996        1.22    0.82
#> areaFM:year_f1996        2.76    0.86
#> areaHO:year_f1996        0.12    0.89
#> areaJM:year_f1996        0.67    0.88
#> areaMU:year_f1996       -1.12    0.81
#> areaRA:year_f1996        0.08    0.91
#> areaSI_EK:year_f1996    -1.63    0.82
#> areaSI_HA:year_f1996     4.98    0.82
#> areaTH:year_f1996        0.67    0.93
#> areaVN:year_f1996       -1.91    0.83
#> areaBT:year_f1997        4.81    0.81
#> areaFB:year_f1997       -1.48    0.83
#> areaFM:year_f1997       -1.75    0.86
#> areaHO:year_f1997       -2.61    0.96
#> areaJM:year_f1997       -0.56    0.88
#> areaMU:year_f1997       -3.61    0.80
#> areaRA:year_f1997       -1.83    0.92
#> areaSI_EK:year_f1997    -2.86    0.82
#> areaSI_HA:year_f1997     3.38    0.83
#> areaTH:year_f1997       -1.91    0.92
#> areaVN:year_f1997        0.54    0.78
#> areaBT:year_f1998        8.71    0.79
#> areaFB:year_f1998        1.10    0.79
#> areaFM:year_f1998        1.38    0.86
#> areaHO:year_f1998       -0.18    0.87
#> areaJM:year_f1998        1.16    0.82
#> areaMU:year_f1998       -1.04    0.78
#> areaRA:year_f1998        0.17    0.86
#> areaSI_EK:year_f1998    -1.32    0.81
#> areaSI_HA:year_f1998     6.26    0.79
#> areaTH:year_f1998        1.47    0.89
#> areaVN:year_f1998       -5.64    0.77
#> areaBT:year_f1999        7.94    0.79
#> areaFB:year_f1999        1.32    0.80
#> areaFM:year_f1999        0.69    0.85
#> areaHO:year_f1999       -2.95    0.93
#> areaJM:year_f1999        1.16    0.82
#> areaMU:year_f1999       -1.55    0.78
#> areaRA:year_f1999       -1.43    0.87
#> areaSI_EK:year_f1999    -2.25    0.80
#> areaSI_HA:year_f1999     5.00    0.79
#> areaTH:year_f1999        0.80    0.89
#> areaVN:year_f1999       -1.78    0.78
#> areaBT:year_f2000        6.76    0.79
#> areaFB:year_f2000        1.28    0.80
#> areaFM:year_f2000        1.08    0.84
#> areaHO:year_f2000       -0.26    0.87
#> areaJM:year_f2000        1.80    0.82
#> areaMU:year_f2000       -0.52    0.77
#> areaRA:year_f2000       -0.60    0.86
#> areaSI_EK:year_f2000    -0.71    0.79
#> areaSI_HA:year_f2000     5.89    0.79
#> areaTH:year_f2000        0.90    0.89
#> areaVN:year_f2000        0.44    0.76
#> areaBT:year_f2001        7.54    0.79
#> areaFB:year_f2001        0.79    0.80
#> areaFM:year_f2001       -0.02    0.84
#> areaHO:year_f2001       -1.29    0.87
#> areaJM:year_f2001        0.53    0.82
#> areaMU:year_f2001       -2.12    0.77
#> areaRA:year_f2001       -0.38    0.86
#> areaSI_EK:year_f2001    -2.90    0.80
#> areaSI_HA:year_f2001     4.65    0.79
#> areaTH:year_f2001        0.11    0.89
#> areaVN:year_f2001       -4.57    0.78
#> areaBT:year_f2002        5.41    0.79
#> areaFB:year_f2002        0.70    0.81
#> areaFM:year_f2002        1.64    0.84
#> areaHO:year_f2002       -1.85    0.87
#> areaJM:year_f2002        0.72    0.82
#> areaMU:year_f2002       -2.99    0.78
#> areaRA:year_f2002       -0.81    0.87
#> areaSI_EK:year_f2002    -1.90    0.80
#> areaSI_HA:year_f2002     4.16    0.79
#> areaTH:year_f2002        0.31    0.83
#> areaVN:year_f2002        0.44    0.76
#> areaBT:year_f2003        5.36    0.80
#> areaFB:year_f2003        0.11    0.81
#> areaFM:year_f2003        0.73    0.85
#> areaHO:year_f2003       -0.77    0.88
#> areaJM:year_f2003       -0.17    0.83
#> areaMU:year_f2003       -2.91    0.78
#> areaRA:year_f2003       -1.20    0.88
#> areaSI_EK:year_f2003    -2.61    0.81
#> areaSI_HA:year_f2003     2.92    0.80
#> areaTH:year_f2003       -1.10    0.84
#> areaVN:year_f2003        1.27    0.77
#> areaBT:year_f2004        6.47    0.79
#> areaFB:year_f2004        1.00    0.80
#> areaFM:year_f2004        0.94    0.84
#> areaHO:year_f2004       -0.29    0.87
#> areaJM:year_f2004        1.04    0.82
#> areaMU:year_f2004       -1.63    0.77
#> areaRA:year_f2004       -0.90    0.87
#> areaSI_EK:year_f2004    -1.36    0.80
#> areaSI_HA:year_f2004     5.52    0.79
#> areaTH:year_f2004        1.14    0.82
#> areaVN:year_f2004        2.38    0.77
#> areaBT:year_f2005        7.07    0.79
#> areaFB:year_f2005        0.71    0.80
#> areaFM:year_f2005        0.17    0.84
#> areaHO:year_f2005       -1.08    0.87
#> areaJM:year_f2005        0.41    0.82
#> areaMU:year_f2005       -2.21    0.78
#> areaRA:year_f2005       -0.56    0.87
#> areaSI_EK:year_f2005    -2.36    0.83
#> areaSI_HA:year_f2005     4.39    0.79
#> areaTH:year_f2005       -0.24    0.83
#> areaVN:year_f2005       -1.37    0.76
#> areaBT:year_f2006        7.61    0.92
#> areaFB:year_f2006        3.36    0.84
#> areaFM:year_f2006        1.33    0.87
#> areaHO:year_f2006        1.74    0.90
#> areaJM:year_f2006        3.87    0.85
#> areaMU:year_f2006        0.27    0.81
#> areaRA:year_f2006        3.58    1.00
#> areaSI_EK:year_f2006     1.47    0.83
#> areaSI_HA:year_f2006     7.10    0.83
#> areaTH:year_f2006        0.87    0.78
#> areaVN:year_f2006        2.58    0.80
#> areaBT:year_f2007        8.89    0.84
#> areaFB:year_f2007        1.68    0.83
#> areaFM:year_f2007        0.97    0.87
#> areaHO:year_f2007        0.24    0.90
#> areaJM:year_f2007        1.78    0.85
#> areaMU:year_f2007       -0.39    0.81
#> areaRA:year_f2007        0.16    0.90
#> areaSI_EK:year_f2007    -0.30    0.83
#> areaSI_HA:year_f2007     6.78    0.82
#> areaTH:year_f2007        1.50    0.84
#> areaVN:year_f2007        1.39    0.83
#> areaBT:year_f2008        7.80    0.87
#> areaFB:year_f2008        1.09    0.83
#> areaFM:year_f2008        0.84    0.87
#> areaHO:year_f2008       -0.92    0.90
#> areaJM:year_f2008        1.09    0.85
#> areaMU:year_f2008       -0.65    0.81
#> areaRA:year_f2008       -0.96    0.90
#> areaSI_EK:year_f2008    -0.51    0.83
#> areaSI_HA:year_f2008     6.41    0.82
#> areaTH:year_f2008       -0.24    0.87
#> areaVN:year_f2008       -0.41    0.82
#> areaBT:year_f2009        5.62    0.90
#> areaFB:year_f2009        0.73    0.83
#> areaFM:year_f2009       -0.12    0.87
#> areaHO:year_f2009       -0.92    0.90
#> areaJM:year_f2009        0.03    0.84
#> areaMU:year_f2009       -1.74    0.80
#> areaRA:year_f2009        0.25    0.90
#> areaSI_EK:year_f2009    -2.69    0.83
#> areaSI_HA:year_f2009     3.84    0.82
#> areaTH:year_f2009       -0.27    0.84
#> areaVN:year_f2009        0.79    0.84
#> areaBT:year_f2010        3.51    0.99
#> areaFB:year_f2010        0.48    0.95
#> areaFM:year_f2010        0.00    0.99
#> areaHO:year_f2010       -0.95    1.01
#> areaJM:year_f2010        0.72    0.97
#> areaMU:year_f2010       -1.61    0.93
#> areaRA:year_f2010        0.18    1.01
#> areaSI_EK:year_f2010    -0.33    1.02
#> areaSI_HA:year_f2010     4.22    0.95
#> areaTH:year_f2010       -3.16    0.92
#> areaVN:year_f2010        1.66    0.98
#> areaBT:year_f2011        5.94    1.02
#> areaFB:year_f2011        1.28    0.95
#> areaFM:year_f2011        0.80    0.99
#> areaHO:year_f2011        0.14    1.01
#> areaJM:year_f2011        0.80    0.97
#> areaMU:year_f2011       -1.30    0.93
#> areaRA:year_f2011        0.64    1.01
#> areaSI_EK:year_f2011    -1.61    0.95
#> areaSI_HA:year_f2011     4.16    0.95
#> areaTH:year_f2011       -2.72    0.90
#> areaVN:year_f2011        0.28    0.96
#> areaBT:year_f2012        6.16    0.98
#> areaFB:year_f2012        0.39    0.93
#> areaFM:year_f2012        0.26    0.96
#> areaHO:year_f2012       -1.06    0.99
#> areaJM:year_f2012       -0.02    0.94
#> areaMU:year_f2012       -1.78    0.91
#> areaRA:year_f2012       -0.36    1.04
#> areaSI_EK:year_f2012    -2.13    0.92
#> areaSI_HA:year_f2012     3.29    0.92
#> areaTH:year_f2012       -1.97    0.87
#> areaVN:year_f2012       -2.88    0.90
#> areaBT:year_f2013        5.14    1.00
#> areaFB:year_f2013        0.94    0.94
#> areaFM:year_f2013        0.75    0.98
#> areaHO:year_f2013       -0.33    1.01
#> areaJM:year_f2013        0.84    0.96
#> areaMU:year_f2013       -1.18    0.92
#> areaRA:year_f2013       -0.74    1.05
#> areaSI_EK:year_f2013    -1.94    0.94
#> areaSI_HA:year_f2013     1.84    0.94
#> areaTH:year_f2013       -1.27    0.89
#> areaVN:year_f2013        0.91    0.91
#> areaBT:year_f2014        5.07    0.97
#> areaFB:year_f2014        0.82    0.94
#> areaFM:year_f2014        0.12    0.97
#> areaHO:year_f2014       -0.74    1.00
#> areaJM:year_f2014        1.27    0.95
#> areaMU:year_f2014       -0.85    0.92
#> areaRA:year_f2014        0.42    1.00
#> areaSI_EK:year_f2014    -1.21    0.93
#> areaSI_HA:year_f2014     3.08    0.93
#> areaTH:year_f2014       -0.84    0.88
#> areaVN:year_f2014        1.32    0.96
#> areaBT:year_f2015        6.48    1.00
#> areaFB:year_f2015        0.09    0.93
#> areaFM:year_f2015       -0.34    0.96
#> areaHO:year_f2015       -1.67    1.00
#> areaJM:year_f2015        0.41    0.94
#> areaMU:year_f2015       -1.29    0.91
#> areaRA:year_f2015       -1.19    0.99
#> areaSI_EK:year_f2015    -1.69    0.92
#> areaSI_HA:year_f2015     2.22    0.92
#> areaTH:year_f2015       -1.10    0.88
#> areaVN:year_f2015       -0.23    0.94
#> areaBT:year_f2016        5.40    0.95
#> areaFB:year_f2016        0.19    0.94
#> areaFM:year_f2016       -1.16    0.97
#> areaHO:year_f2016       -1.75    0.99
#> areaJM:year_f2016       -0.03    0.95
#> areaMU:year_f2016       -2.01    0.91
#> areaRA:year_f2016       -0.54    0.99
#> areaSI_EK:year_f2016    -2.18    0.93
#> areaSI_HA:year_f2016     2.92    0.93
#> areaTH:year_f2016       -1.50    0.88
#> areaVN:year_f2016       -2.88    1.26
#> areaBT:year_f2017        5.82    0.94
#> areaFB:year_f2017       -0.11    0.93
#> areaFM:year_f2017       -0.70    0.96
#> areaHO:year_f2017       -1.73    0.98
#> areaJM:year_f2017        0.21    0.94
#> areaMU:year_f2017       -1.60    0.91
#> areaRA:year_f2017       -1.44    0.99
#> areaSI_EK:year_f2017    -1.84    0.92
#> areaSI_HA:year_f2017     2.55    0.92
#> areaTH:year_f2017       -2.08    0.87
#> areaVN:year_f2017       -1.17    0.95
#> areaBT:year_f2018        5.04    0.96
#> areaFB:year_f2018        0.60    0.94
#> areaFM:year_f2018        0.73    0.97
#> areaHO:year_f2018       -1.86    0.99
#> areaJM:year_f2018        0.53    0.96
#> areaMU:year_f2018       -2.43    0.92
#> areaRA:year_f2018       -0.86    1.00
#> areaSI_EK:year_f2018    -2.56    0.93
#> areaSI_HA:year_f2018     2.73    0.93
#> areaTH:year_f2018       -0.52    0.94
#> areaVN:year_f2018       -0.92    0.97
#> areaBT:year_f2019        5.74    1.03
#> areaFB:year_f2019       -0.50    1.01
#> areaFM:year_f2019        0.00    0.96
#> areaHO:year_f2019       -1.43    0.99
#> areaJM:year_f2019        0.45    0.95
#> areaMU:year_f2019       -1.81    0.91
#> areaRA:year_f2019       -0.96    0.99
#> areaSI_EK:year_f2019    -1.34    0.92
#> areaSI_HA:year_f2019     4.14    0.92
#> areaTH:year_f2019       -1.56    0.88
#> areaVN:year_f2019        0.10    0.99
#> areaBT:year_f2020        7.64    0.96
#> areaFB:year_f2020       -0.71    1.00
#> areaFM:year_f2020       -1.48    0.98
#> areaHO:year_f2020       -1.11    1.15
#> areaJM:year_f2020        0.25    0.98
#> areaMU:year_f2020       -0.60    0.93
#> areaRA:year_f2020       -0.81    1.05
#> areaSI_EK:year_f2020    -1.91    1.13
#> areaSI_HA:year_f2020     2.19    0.99
#> areaTH:year_f2020        0.92    0.95
#> areaVN:year_f2020       -0.49    0.93
#> areaBT:year_f2021        4.45    1.12
#> areaFB:year_f2021        0.26    1.08
#> areaFM:year_f2021        0.12    1.03
#> areaHO:year_f2021       -1.25    1.12
#> areaJM:year_f2021        0.57    0.98
#> areaMU:year_f2021        0.66    0.95
#> areaRA:year_f2021       -2.40    1.07
#> areaSI_EK:year_f2021    -0.56    1.10
#> areaSI_HA:year_f2021     2.03    1.01
#> areaTH:year_f2021       -1.29    0.99
#> areaVN:year_f2021       -0.11    0.99
#> areaBT:year_f2022        7.22    1.51
#> areaFB:year_f2022        0.29    1.84
#> areaFM:year_f2022       -0.11    1.50
#> areaHO:year_f2022       -1.60    1.60
#> areaJM:year_f2022       -0.28    1.51
#> areaMU:year_f2022        0.78    1.42
#> areaRA:year_f2022       -1.72    1.53
#> areaSI_EK:year_f2022     0.15    1.68
#> areaSI_HA:year_f2022     3.24    1.48
#> areaTH:year_f2022        0.66    1.47
#> areaVN:year_f2022       -4.72    1.46
#> 
#> Smooth terms:
#>                      Std. Dev.
#> sds(yday):areaBS)         0.89
#> sds(yday):areaBT)         0.66
#> sds(yday):areaFB)         1.04
#> sds(yday):areaFM)         0.96
#> sds(yday):areaHO)         1.16
#> sds(yday):areaJM)         1.34
#> sds(yday):areaMU)         0.88
#> sds(yday):areaRA)         1.40
#> sds(yday):areaSI_EK)      0.73
#> sds(yday):areaSI_HA)      0.73
#> sds(yday):areaTH)         0.78
#> sds(yday):areaVN)         0.78
#> 
#> Dispersion parameter: 1.80
#> ML criterion at convergence: 219571.168
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

unique(d$area)
#>  [1] BS    BT    FB    FM    HO    JM    MU    RA    SI_EK SI_HA TH    VN   
#> Levels: BS BT FB FM HO JM MU RA SI_EK SI_HA TH VN
unique(d$source_f)
#> [1] logger  errs    fishing
#> Levels: errs fishing logger

Predict

Code
# Make a new data frame and predict!
nd <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                             area = unique(d$area),
                             source = unique(d$source_f),
                             year = unique(d$year))) |>
  mutate(source_f = as.factor(source),
         year_f = as.factor(year)) |> # Left join in growth data column
  left_join(gdat, by = "area") |> 
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x           0
#>            > rows only in y  (        0)
#>            > matched rows     1,093,608
#>            >                 ===========
#>            > rows total       1,093,608
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA

# Predict
nd$pred <- predict(m, newdata = nd)$est

In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear

Code
nd <- nd |> 
  mutate(keep = "Y",
         keep = ifelse(area == "BT" & year < 1980, "N", keep),
         keep = ifelse(area == "SI_HA" & year < 1972, "N", keep)) |>
  filter(keep == "Y")
#> mutate: new variable 'keep' (character) with 2 unique values and 0% NA
#> filter: removed 79,056 rows (7%), 1,014,552 rows remaining
  
nd_sub <- nd |> 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year < 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year < 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y")
#> mutate: changed 935,496 values (92%) of 'keep' (0 new NA)
#> filter: removed 935,496 rows (92%), 79,056 rows remaining

# Now change the labels to BT and SI_EK...
nd_sub <- nd_sub |> 
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
#> mutate: converted 'area' from factor to character (0 new NA)

# Bind rows and plot only the temperature series we will use for growth modelling
nd <- bind_rows(nd, nd_sub) |> select(-keep)
#> select: dropped one variable (keep)

Plot predictions

Code
# Trimmed plot
nd |> 
  filter(growth_dat == "Y") |> 
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
    geom_raster() +
    facet_wrap(~area, ncol = 3) +
    scale_fill_viridis(option = "magma") +
    scale_color_viridis(option = "magma") +
    labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)")
#> filter: removed 648,918 rows (59%), 444,690 rows remaining

Code

# Full plot
nd |>
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
    geom_raster() +
    facet_wrap(~area, ncol = 3) +
    scale_fill_viridis(option = "magma") +
    scale_color_viridis(option = "magma") +
    labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)")

Detailed exploration of predictions

Code
# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

for(i in unique(nd$area)) {
  
  plotdat <- nd |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = c(0.8, 0.08), 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/full_model/temp_pred_yday_area_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 95,167 rows (97%), 3,449 rows remaining
#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 89,164 rows (90%), 9,452 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 89,155 rows (90%), 9,461 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 86,660 rows (88%), 11,956 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 90,466 rows (92%), 8,150 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 88,494 rows (90%), 10,122 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 87,743 rows (89%), 10,873 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 93,808 rows (95%), 4,808 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 89,755 rows (91%), 8,861 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 82,797 rows (84%), 15,819 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 95,212 rows (97%), 3,404 rows remaining

#> filter: removed 1,002,474 rows (92%), 91,134 rows remaining
#> filter: removed 96,499 rows (98%), 2,117 rows remaining

Area-specific models

Code
spec_preds <- list()

# Drop VN, no logger data? 
dspec <- d |> filter(!area == "VN")
#> filter: removed 2,129 rows (2%), 96,487 rows remaining

for(i in unique(dspec$area)) {
  
  dd <- dspec |> filter(area == i)

  mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                 data = dd,
                 family = student(df = 8),
                 spatial = "off",
                 spatiotemporal = "off",
                 knots = list(yday = c(0.5, 364.5)),
                 control = sdmTMBcontrol(newton_loops = 1))
  
  sanity(mspec)

  # QQ plots
  mcmc_res_msep <- residuals(mspec, type = "mle-mcmc",
                             mcmc_samples = sdmTMBextra::predict_mle_mcmc(mspec,
                                                                          mcmc_iter = 201,
                                                                          mcmc_warmup = 200))
  
  print(ggplot(dd, aes(sample = mcmc_res_msep)) +
    stat_qq(size = 0.75, shape = 21, fill = NA) +
    stat_qq_line() +
    labs(y = "Sample Quantiles", x = "Theoretical Quantiles", title = paste("Area = ", i)) + 
    theme(aspect.ratio = 1))
  
  ggsave(paste(home, "/figures/supp/qq_temp", i, ".pdf", sep = ""), width = 11, height = 11, units = "cm")

  # Predict on new data
  nd_area <- data.frame(expand.grid(yday = seq(min(dd$yday), max(dd$yday), by = 1),
                                    area = unique(dd$area),
                                    source = unique(dd$source_f),
                                    year = unique(dd$year))) |>
    mutate(source_f = as.factor(source),
           year_f = as.factor(year)) |> 
    left_join(gdat, by = "area") |> 
    mutate(area = as.factor(area),
           growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
  
  nd_area$pred <- predict(mspec, newdata = nd_area)$est

  # Save!
  spec_preds[[i]] <- nd_area
  
}
#> filter: removed 93,026 rows (96%), 3,461 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000694 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.94 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.22539 seconds (Warm-up)
#> Chain 1:                0.003849 seconds (Sampling)
#> Chain 1:                3.22924 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     83,664
#>            >                 ========
#>            > rows total       83,664
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 87,023 rows (90%), 9,464 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001521 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 15.21 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 7.44321 seconds (Warm-up)
#> Chain 1:                0.010882 seconds (Sampling)
#> Chain 1:                7.45409 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     91,134
#>            >                 ========
#>            > rows total       91,134
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 87,014 rows (90%), 9,473 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001426 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.26 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 7.65687 seconds (Warm-up)
#> Chain 1:                0.021189 seconds (Sampling)
#> Chain 1:                7.67805 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     90,885
#>            >                 ========
#>            > rows total       90,885
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 84,519 rows (88%), 11,968 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001655 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16.55 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 8.30525 seconds (Warm-up)
#> Chain 1:                0.012321 seconds (Sampling)
#> Chain 1:                8.31757 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     91,134
#>            >                 ========
#>            > rows total       91,134
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 88,325 rows (92%), 8,162 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001278 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.78 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 7.69745 seconds (Warm-up)
#> Chain 1:                0.018345 seconds (Sampling)
#> Chain 1:                7.71579 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     90,885
#>            >                 ========
#>            > rows total       90,885
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 86,353 rows (89%), 10,134 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.00153 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 15.3 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 7.65834 seconds (Warm-up)
#> Chain 1:                0.011259 seconds (Sampling)
#> Chain 1:                7.6696 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     84,660
#>            >                 ========
#>            > rows total       84,660
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 85,602 rows (89%), 10,885 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001713 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 17.13 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 7.66399 seconds (Warm-up)
#> Chain 1:                0.012704 seconds (Sampling)
#> Chain 1:                7.67669 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     91,134
#>            >                 ========
#>            > rows total       91,134
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 91,667 rows (95%), 4,820 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000763 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.63 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.80462 seconds (Warm-up)
#> Chain 1:                0.011394 seconds (Sampling)
#> Chain 1:                3.81601 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     91,134
#>            >                 ========
#>            > rows total       91,134
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 87,614 rows (91%), 8,873 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001443 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.43 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 6.26786 seconds (Warm-up)
#> Chain 1:                0.010224 seconds (Sampling)
#> Chain 1:                6.27809 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     91,134
#>            >                 ========
#>            > rows total       91,134
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 80,656 rows (84%), 15,831 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.002412 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 24.12 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 9.63288 seconds (Warm-up)
#> Chain 1:                0.017792 seconds (Sampling)
#> Chain 1:                9.65068 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     91,134
#>            >                 ========
#>            > rows total       91,134
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA
#> filter: removed 93,071 rows (96%), 3,416 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf

#> Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
#> Inf
#> Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
#> -Inf

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000584 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.84 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 2.85052 seconds (Warm-up)
#> Chain 1:                0.004189 seconds (Sampling)
#> Chain 1:                2.85471 seconds (Total)
#> Chain 1:
#> mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
#>         new variable 'year_f' (factor) with 83 unique values and 0% NA
#> left_join: added 2 columns (min, max)
#>            > rows only in x        0
#>            > rows only in y  (    11)
#>            > matched rows     90,885
#>            >                 ========
#>            > rows total       90,885
#> mutate: converted 'area' from character to factor (0 new NA)
#>         new variable 'growth_dat' (character) with 2 unique values and 0% NA

Code

nd_area <- dplyr::bind_rows(spec_preds)

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear

nd_area <- nd_area |> 
  mutate(keep = "Y",
         keep = ifelse(area == "BT" & year < 1980, "N", keep),
         keep = ifelse(area == "SI_HA" & year < 1972, "N", keep)) |>
  filter(keep == "Y")
#> mutate: new variable 'keep' (character) with 2 unique values and 0% NA
#> filter: removed 79,056 rows (8%), 908,727 rows remaining
  
nd_area_sub <- nd_area |> 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year < 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year < 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y")
#> mutate: changed 829,671 values (91%) of 'keep' (0 new NA)
#> filter: removed 829,671 rows (91%), 79,056 rows remaining

# Now change the labels to BT and SI_EK...
nd_area_sub <- nd_area_sub |> 
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
#> mutate: converted 'area' from factor to character (0 new NA)

# Bind rows and plot only the temperature series we will use for growth modelling
nd_area <- bind_rows(nd_area, nd_area_sub) |> select(-keep)
#> select: dropped one variable (keep)

nd_area |> 
  filter(area %in% c("FM", "BT", "SI_EK", "SI_HA")) |> 
  filter(year <= 1980 & year >= 1966) |> 
  group_by(area, year) |> 
  summarise(mean_temp = mean(pred)) |> 
  ungroup() |> 
  pivot_wider(names_from = area, values_from = mean_temp)
#> filter: removed 623,247 rows (63%), 364,536 rows remaining
#> filter: removed 298,656 rows (82%), 65,880 rows remaining
#> group_by: 2 grouping variables (area, year)
#> summarise: now 60 rows and 3 columns, one group variable remaining (area)
#> ungroup: no grouping variables
#> pivot_wider: reorganized (area, mean_temp) into (BT, FM, SI_EK, SI_HA) [was 60x3, now 15x5]
#> # A tibble: 15 × 5
#>     year    BT    FM SI_EK SI_HA
#>    <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1  1966  8.23  8.23  7.96  7.96
#>  2  1967  9.05  9.05  8.74  8.74
#>  3  1968  8.74  8.74  8.36  8.36
#>  4  1969  8.60  8.60  8.48  8.48
#>  5  1970  8.38  8.38  8.02  8.02
#>  6  1971  8.40  8.40  8.32  8.32
#>  7  1972  8.86  8.86  8.87  9.82
#>  8  1973  9.24  9.24  9.00 11.2 
#>  9  1974  9.13  9.13  8.65  8.43
#> 10  1975  7.40  7.40  9.30 14.6 
#> 11  1976  7.01  7.01  7.97 13.7 
#> 12  1977  5.60  5.60  7.79 13.0 
#> 13  1978  5.55  5.55  7.57 13.1 
#> 14  1979  5.82  5.82  7.40 11.3 
#> 15  1980 11.9   6.87  7.71 13.4

Plot detailed exploration of predictions

Code
# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

for(i in unique(nd_area$area)) {
  
  plotdat <- nd_area |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
#> filter: removed 904,119 rows (92%), 83,664 rows remaining
#> filter: removed 95,167 rows (97%), 3,449 rows remaining
#> filter: removed 896,649 rows (91%), 91,134 rows remaining
#> filter: removed 89,164 rows (90%), 9,452 rows remaining

#> filter: removed 896,898 rows (91%), 90,885 rows remaining
#> filter: removed 89,155 rows (90%), 9,461 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining
#> filter: removed 86,660 rows (88%), 11,956 rows remaining

#> filter: removed 896,898 rows (91%), 90,885 rows remaining
#> filter: removed 90,466 rows (92%), 8,150 rows remaining

#> filter: removed 903,123 rows (91%), 84,660 rows remaining
#> filter: removed 88,494 rows (90%), 10,122 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining
#> filter: removed 87,743 rows (89%), 10,873 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining
#> filter: removed 93,808 rows (95%), 4,808 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining
#> filter: removed 89,755 rows (91%), 8,861 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining
#> filter: removed 82,797 rows (84%), 15,819 rows remaining

#> filter: removed 896,898 rows (91%), 90,885 rows remaining
#> filter: removed 95,212 rows (97%), 3,404 rows remaining

Compare predictions from full and area-specific models

Code
long_preds <- bind_rows(nd_area |> filter(growth_dat == "Y") |> mutate(model = "area-specific"),
                        nd |> filter(growth_dat == "Y") |> mutate(model = "full model"))
#> filter: removed 563,175 rows (57%), 424,608 rows remaining
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 648,918 rows (59%), 444,690 rows remaining
#> mutate: new variable 'model' (character) with one unique value and 0% NA


for(i in unique(long_preds$area)) {
  
  plotdat <- long_preds |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source, linetype = model)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & growth_dat == "Y"), size = 0.2,
                 aes(yday, temp, color = source), inherit.aes = FALSE) + 
      geom_line(linewidth = 0.3) + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)",
           title = paste("Area = ", i), color = "", linetype = "")
  )
  
  ggsave(paste0(home, "/figures/supp/full_model/temp_pred_comp_area_model", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
#> filter: removed 831,390 rows (96%), 37,908 rows remaining
#> filter: removed 96,724 rows (98%), 1,892 rows remaining
#> filter: removed 785,850 rows (90%), 83,448 rows remaining
#> filter: removed 92,124 rows (93%), 6,492 rows remaining

#> filter: removed 764,034 rows (88%), 105,264 rows remaining
#> filter: removed 90,040 rows (91%), 8,576 rows remaining

#> filter: removed 785,850 rows (90%), 83,448 rows remaining
#> filter: removed 92,321 rows (94%), 6,295 rows remaining

#> filter: removed 792,543 rows (91%), 76,755 rows remaining
#> filter: removed 91,764 rows (93%), 6,852 rows remaining

#> filter: removed 733,746 rows (84%), 135,552 rows remaining
#> filter: removed 89,689 rows (91%), 8,927 rows remaining

#> filter: removed 825,378 rows (95%), 43,920 rows remaining
#> filter: removed 96,656 rows (98%), 1,960 rows remaining

#> filter: removed 820,986 rows (94%), 48,312 rows remaining
#> filter: removed 96,380 rows (98%), 2,236 rows remaining

#> filter: removed 763,890 rows (88%), 105,408 rows remaining
#> filter: removed 91,638 rows (93%), 6,978 rows remaining

#> filter: removed 766,086 rows (88%), 103,212 rows remaining
#> filter: removed 85,170 rows (86%), 13,446 rows remaining

#> filter: removed 836,403 rows (96%), 32,895 rows remaining
#> filter: removed 96,832 rows (98%), 1,784 rows remaining

#> filter: removed 856,122 rows (98%), 13,176 rows remaining
#> filter: removed 98,238 rows (>99%), 378 rows remaining

Code


# Join wide data to plot differences easily
wide_pred <- left_join(nd_area |> 
                         dplyr::select(yday, area, year, pred, source) |> 
                         pivot_wider(names_from = source, values_from = pred) |> 
                         rename(logger_as = logger,
                                errs_as = errs,
                                fishing_as = fishing),
                       
                       nd |> 
                         dplyr::select(yday, area, year, pred, source) |> 
                         pivot_wider(names_from = source, values_from = pred) |> 
                         rename(logger_full = logger,
                                errs_full = errs,
                                fishing_full = fishing),
                       
                       by = c("yday", "area", "year"))
#> pivot_wider: reorganized (pred, source) into (logger, errs, fishing) [was 987783x5, now 329261x6]
#> rename: renamed 3 variables (logger_as, errs_as, fishing_as)
#> pivot_wider: reorganized (pred, source) into (logger, errs, fishing) [was 1093608x5, now 364536x6]
#> rename: renamed 3 variables (logger_full, errs_full, fishing_full)
#> left_join: added 3 columns (logger_full, errs_full, fishing_full)
#>            > rows only in x         0
#>            > rows only in y  ( 35,275)
#>            > matched rows     329,261
#>            >                 =========
#>            > rows total       329,261
  
  
diff_pred <- wide_pred |> 
  mutate(logger_diff = logger_as - logger_full,
         errs_diff = errs_as - errs_full,
         fishing_diff = fishing_as - fishing_full) |> 
  select(yday, area, year, logger_diff, errs_diff, fishing_diff) |> 
  pivot_longer(c("logger_diff", "errs_diff", "fishing_diff"))
#> mutate: new variable 'logger_diff' (double) with 302,909 unique values and 0% NA
#>         new variable 'errs_diff' (double) with 302,909 unique values and 0% NA
#>         new variable 'fishing_diff' (double) with 302,909 unique values and 0% NA
#> select: dropped 6 variables (logger_as, errs_as, fishing_as, logger_full, errs_full, …)
#> pivot_longer: reorganized (logger_diff, errs_diff, fishing_diff) into (name, value) [was 329261x6, now 987783x5]


for(i in unique(diff_pred$area)) {
  
  plotdat <- diff_pred |> filter(area == i)
  
  print(ggplot(diff_pred, aes(yday, value, color = name)) + 
      facet_wrap(~year) +  
      scale_color_brewer(palette = "Dark2") + 
      geom_hline(yintercept = 0, linetype = 2, alpha = 0.5) + 
      geom_line() +
      labs(x = "Day of the year", y = "Predicted SST (°C)",
           title = paste("Area = ", i), color = "", linetype = ""))
}
#> filter: removed 904,119 rows (92%), 83,664 rows remaining
#> filter: removed 896,649 rows (91%), 91,134 rows remaining

#> filter: removed 896,898 rows (91%), 90,885 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining

#> filter: removed 896,898 rows (91%), 90,885 rows remaining

#> filter: removed 903,123 rows (91%), 84,660 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining

#> filter: removed 896,649 rows (91%), 91,134 rows remaining

#> filter: removed 896,898 rows (91%), 90,885 rows remaining

Plot summarized data and predictions

Code
# Summarise data
dsum <- d |> 
  group_by(year, area, source) |> 
  summarise(temp = mean(temp)) |> 
  mutate(type = "data")
#> group_by: 3 grouping variables (year, area, source)
#> summarise: now 1,676 rows and 4 columns, 2 group variables remaining (year, area)
#> mutate (grouped): new variable 'type' (character) with one unique value and 0% NA

# Summarise predictions from full model and area-specific model
preds <- nd |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  group_by(area, year) |> 
  summarise(temp = mean(pred)) |> 
  mutate(model = "full model")
#> filter: removed 945,378 rows (86%), 148,230 rows remaining
#> group_by: 2 grouping variables (area, year)
#> summarise: now 405 rows and 3 columns, one group variable remaining (area)
#> mutate (grouped): new variable 'model' (character) with one unique value and 0% NA

preds_area <- nd_area |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  group_by(area, year) |> 
  summarise(temp = mean(pred)) |> 
  mutate(model = "area model")
#> filter: removed 846,247 rows (86%), 141,536 rows remaining
#> group_by: 2 grouping variables (area, year)
#> summarise: now 393 rows and 3 columns, one group variable remaining (area)
#> mutate (grouped): new variable 'model' (character) with one unique value and 0% NA

preds_comb <- bind_rows(preds, preds_area)

ggplot(preds_comb, aes(year, temp, color = source, linetype = model)) + 
  geom_point(data = dsum, aes(year, temp, color = source), size = 0.75, alpha = 0.75, inherit.aes = FALSE) + 
  scale_color_brewer(palette = "Accent") +
  geom_line(linewidth = 0.5, color = "grey20") + 
  facet_wrap(~area) +
  theme(legend.position = "bottom")

Make final plot using the area-specific model

Code
order <- preds_area |>
  group_by(area) |> 
  summarise(mean_temp = mean(temp)) |> 
  arrange(desc(mean_temp))
#> group_by: one grouping variable (area)
#> summarise: now 11 rows and 2 columns, ungrouped

order
#> # A tibble: 11 × 2
#>    area  mean_temp
#>    <chr>     <dbl>
#>  1 SI_HA     14.1 
#>  2 BT        12.2 
#>  3 TH         9.30
#>  4 FM         8.93
#>  5 JM         8.84
#>  6 BS         8.55
#>  7 SI_EK      8.27
#>  8 FB         8.20
#>  9 MU         6.96
#> 10 RA         6.61
#> 11 HO         6.60
# Save plot order..

topt <- 11.8 # Overall t_opt from 02-fit-vbge.qmd! Update if needed

# Add latitude
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat","lon"), as.numeric) |> 
  arrange(desc(lat))
#> mutate_at: converted 'lat' from character to double (0 new NA)
#>            converted 'lon' from character to double (0 new NA)

ggplot(preds_area, aes(year, temp, color = temp)) + 
  facet_wrap(~factor(area, levels = area_attr$area), ncol = 3) + 
  geom_hline(yintercept = topt, linewidth = 0.3, linetype = 2, color = "grey") +
  geom_line() +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_viridis(option = "magma", name = "Area") +
  guides(color = "none") 

Code

preds |> 
  group_by(area) |> 
  summarise(min = min(year),
            max = max(year)) |> 
  arrange(min)
#> group_by: one grouping variable (area)
#> summarise: now 12 rows and 3 columns, ungrouped
#> # A tibble: 12 × 3
#>    area    min   max
#>    <chr> <dbl> <dbl>
#>  1 JM     1953  2016
#>  2 BT     1963  2000
#>  3 FM     1963  2000
#>  4 SI_EK  1968  2015
#>  5 SI_HA  1968  2014
#>  6 FB     1969  2016
#>  7 MU     1981  2000
#>  8 HO     1982  2016
#>  9 BS     1985  2002
#> 10 RA     1985  2006
#> 11 VN     1987  1998
#> 12 TH     2000  2014

ggsave(paste0(home, "/figures/annual_average_temperature.pdf"), width = 17, height = 17, units = "cm")
Code
# Save prediction df
write_csv(preds, paste0(home, "/output/gam_predicted_temps.csv"))

Code below is not evaluated!

Growing season? This might be different for different areas…

Code
# Find day of the year where temperature exceeds 10C by area across years
# TODO: Also for area-specific predictions?
gs_area <- nd |> 
  group_by(area, yday) |> 
  summarise(mean_pred = mean(pred)) |>
  ungroup() |> 
  filter(mean_pred > 10) |> 
  group_by(area) |> 
  summarise(gs_min = min(yday),
            gs_max = max(yday))
#> group_by: 2 grouping variables (area, yday)
#> summarise: now 4,392 rows and 3 columns, one group variable remaining (area)
#> ungroup: no grouping variables
#> filter: removed 2,556 rows (58%), 1,836 rows remaining
#> group_by: one grouping variable (area)
#> summarise: now 12 rows and 3 columns, ungrouped

nd <- left_join(nd, gs_area, by = "area")
#> left_join: added 2 columns (gs_min, gs_max)
#>            > rows only in x           0
#>            > rows only in y  (        0)
#>            > matched rows     1,093,608
#>            >                 ===========
#>            > rows total       1,093,608

gs_area$mean_pred <- 10

# Plot!
nd |> 
  #filter(growth_dat == "Y") |> 
  group_by(area, yday) |> 
  summarise(mean_pred = mean(pred)) |> 
  ggplot(aes(yday, mean_pred)) +
  geom_line() +
  labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)") + 
  facet_wrap(~area) +
  geom_point(data = gs_area, aes(gs_min, mean_pred), inherit.aes = FALSE, color = "tomato2") + 
  geom_point(data = gs_area, aes(gs_max, mean_pred), inherit.aes = FALSE, color = "tomato2") +
  geom_hline(yintercept = 10, linetype = 2)
#> group_by: 2 grouping variables (area, yday)
#> summarise: now 4,392 rows and 3 columns, one group variable remaining (area)

Now see if there is a systematic pattern in the difference between predicted and observed logger data, which could indicate that the source effect isn’t global but area-specific.

Code
dlog <- d |> 
  filter(source == "logger") |> 
  mutate(type = "data",
         id = paste(area, year, yday, sep = "_")) |> 
  select(id, temp) |> 
  group_by(id) |> 
  summarise(obs = mean(temp)) # sometimes we have more than 1 observation per id
#> filter: removed 34,632 rows (35%), 63,984 rows remaining
#> mutate: changed 63,984 values (100%) of 'id' (63984 fewer NA)
#>         new variable 'type' (character) with one unique value and 0% NA
#> select: dropped 17 variables (year, area, yday, month, date, …)
#> group_by: one grouping variable (id)
#> summarise: now 60,967 rows and 2 columns, ungrouped

# dlog |> 
#   group_by(id) |> 
#   summarise(n = n()) |> 
#   distinct(n)

preds_log <- nd |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  mutate(type = "model",
         id = paste(area, year, yday, sep = "_")) |> 
  filter(id %in% unique(dlog$id)) |> 
  ungroup() |> 
  left_join(dlog, by = "id")
#> filter: removed 945,378 rows (86%), 148,230 rows remaining
#> mutate: new variable 'type' (character) with one unique value and 0% NA
#>         new variable 'id' (character) with 148,230 unique values and 0% NA
#> filter: removed 107,373 rows (72%), 40,857 rows remaining
#> ungroup: no grouping variables
#> left_join: added one column (obs)
#>            > rows only in x        0
#>            > rows only in y  (20,110)
#>            > matched rows     40,857
#>            >                 ========
#>            > rows total       40,857

p1 <- preds_log |> 
  mutate(resid = pred - obs) |> 
  ggplot(aes(as.factor(area), resid, group = as.factor(area))) +
  #geom_jitter(alpha = 0.05, color = "grey20", height = 0, width = 0.2) + 
  geom_violin(fill = "grey70", color = NA) +
  geom_boxplot(width = 0.2, outlier.colour = NA, outlier.color = NA, outlier.fill = NA) +
  guides(color = "none") +
  geom_hline(yintercept = 0, linetype = 2, color = "tomato3", linewidth = 0.75) + 
  labs(x = "Area", y = "Manual residuals")
#> mutate: new variable 'resid' (double) with 40,857 unique values and 0% NA

p1

Code

p1 + facet_wrap(~year) 

Some extra plots for BT especially, because it seems that the offset for logger isn’t that big in years without any logger data.. First predict BT only but with all data sources. Note this code is not run anymore and is instead expanded above to work on all areas

Code
# Predict
ndbt <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                               year = seq(1980, #filter(gdat, area == "BT")$min, 
                                          2004), #filter(gdat, area == "BT")$min),
                               source = unique(d$source))) |>
  mutate(area = "BT") |> 
  mutate(area = as.factor(area),
         source_f = as.factor(source),
         year_f = as.factor(year)) 

ndbt$pred <- predict(m, newdata = ndbt)$est

# Plot
ggplot(ndbt, aes(yday, pred, color = source)) + 
  scale_color_brewer(palette = "Accent") + 
  facet_wrap(~year) + 
  geom_point(data = filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)), size = 0.2,
             aes(yday, temp)) + 
  geom_line() + 
  labs(title = "Area = BT", color = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots.pdf" ), width = 17, height = 17, units = "cm")

# Right, so there is an offset... but it's not big. And the cold temperatures in
# certain years in the 1980's is not due to the offset not working, just that the year
# effect in that year is only informed by cold temperatures and there's no "memory"... 
# Though it seems also that there could be a bigger offset.. perhaps this indicates
# the source should indeed vary by area? Maybe, we if we fit models by area separately instead? Then we don't
# need the area interaction, and we can instead use a random walk or AR1 structure one
# the year effect?

dbt <- d |> filter(area == "BT")

mbt <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
              data = dbt,
              family = student(df = 5),
              spatial = "off",
              spatiotemporal = "off",
              knots = list(yday = c(0.5, 364.5)),
              control = sdmTMBcontrol(newton_loops = 1))

sanity(mbt)
summary(mbt)

# Predict
ndbt2 <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                                year = seq(1980, #filter(gdat, area == "BT")$min, 
                                           2004), #filter(gdat, area == "BT")$min),
                                source = unique(d$source))) |>
  mutate(source_f = as.factor(source),
         year_f = as.factor(year)) 

ndbt2$pred <- predict(mbt, newdata = ndbt2)$est

# Plot without data
ggplot(filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)),
       aes(yday, temp, color = source)) + 
  geom_line(data = ndbt, aes(yday, pred, color = source, linetype = "main model"), alpha = 0.2) +
  geom_line(data = ndbt2, aes(yday, pred, color = source, linetype = "BT specific model"), alpha = 0.2) + 
  scale_color_brewer(palette = "Dark2") + 
  facet_wrap(~year) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots2.pdf" ), width = 17, height = 17, units = "cm")

# Plot with data
ggplot(filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)),
       aes(yday, temp, color = source)) + 
  geom_point(size = 0.001, alpha = 0.3) + 
  geom_line(data = ndbt, aes(yday, pred, color = source, linetype = "main model"), alpha = 0.7) +
  geom_line(data = ndbt2, aes(yday, pred, color = source, linetype = "BT specific model"), alpha = 0.7) + 
  scale_color_brewer(palette = "Dark2") + 
  facet_wrap(~year) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots3.pdf" ), width = 17, height = 17, units = "cm")

# Plot only new predictions, color by year facet by source?
p1 <- ggplot() + 
  geom_line(data = ndbt2, aes(yday, pred, color = year, group = year), alpha = 0.7, linewidth = 0.2) + 
  scale_color_viridis() + 
  facet_wrap(~source) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "Day of the year", y = "Predicted SST (°C)", title = "BT specific model")

p2 <- ggplot() + 
  geom_line(data = ndbt, aes(yday, pred, color = year, group = year), alpha = 0.7, linewidth = 0.2) + 
  scale_color_viridis() + 
  facet_wrap(~source) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "Day of the year", y = "Predicted SST (°C)", title = "main model")

p1 / p2

ggsave(paste0(home, "/figures/supp/BT_test_plots4.pdf" ), width = 17, height = 17, units = "cm")

# Now plot the annual means. Focus on 1980-2004, so that we don't need to worry about
# the Forsmark predictions in years before heating

ndbt_sum <- ndbt |> 
  filter(source == "logger") |> 
  group_by(year) |> 
  summarise(mean_pred = mean(pred)) |> 
  mutate(model = "main model")

ndbt2_sum <- ndbt2 |> 
  filter(source == "logger") |> 
  group_by(year) |> 
  summarise(mean_pred = mean(pred)) |> 
  mutate(model = "BT-specific model")

bt_sums <- bind_rows(ndbt_sum, ndbt2_sum)

ggplot(bt_sums, aes(year, mean_pred, color = model)) + 
  geom_line() +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_brewer(palette = "Accent")

ggsave(paste0(home, "/figures/supp/BT_test_plots5.pdf" ), width = 17, height = 17, units = "cm")

If we want to get uncertainty, we can use nsim instead; this simulates from the linear predictor using the inverse precision matrix, which is a fast way to get a distribution of samples from which we can take e.g. quantiles and means. However, it’s still slow, so the code below isn’t executed yet.

Code
nd_sim <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                                 area = unique(d$area),
                                 year = unique(d$year))) |>
  mutate(source = "logger") |>
  mutate(id = paste(year, area, sep = "_"),
         source_f = as.factor(source),
         year_f = as.factor(year))

# Trim!
nd_sim <- left_join(nd_sim, gdat, by = "area")

nd_sim <- nd_sim |>
  mutate(growth_dat = ifelse(year > min, "Y", "N")) |>
  filter(growth_dat == "Y") |>
  filter(yday %in% c(gs_min:gs_min)) |>
  mutate(area = as.factor(area))

# Predict!
nsim <- 500
m_pred_sims <- predict(m, newdata = nd_sim, nsim = nsim)

# Join sims with prediction data
nd_sim_long <- cbind(nd_sim, as.data.frame(m_pred_sims)) |>
    pivot_longer(c((ncol(nd_sim) + 1):(nsim + ncol(nd_sim))))

# Summarize sims over growing season
sum_pred_gs <- nd_sim_long |>
    ungroup() |>
    group_by(year, area) |>
    summarise(lwr = quantile(value, prob = 0.1),
              est = quantile(value, prob = 0.5),
              upr = quantile(value, prob = 0.9)) |>
    ungroup()

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas..
sum_pred_gs <- preds |>
  mutate(keep = "Y",
         keep = ifelse(area == "BT" & year < 1980, "N", keep),
         keep = ifelse(area == "SI_HA" & year < 1972, "N", keep)) |>
  filter(keep == "Y")

sum_pred_gs_sub <- preds |>
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year < 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year < 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y")

# Now change the labels to BT and SI_EK...
sum_pred_gs_sub <- sum_pred_gs_sub |>
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))

# Bind rows and plot only the temperature series we will use for growth modelling
sum_pred_gs <- bind_rows(sum_pred_gs, sum_pred_gs_sub) |> select(-keep, -type)

order <- sum_pred_gs |>
  group_by(area) |>
  summarise(mean_temp = mean(temp)) |>
  arrange(desc(mean_temp))